%%%%%%%%%%%%
% Calculate investment PLs,
% weekly sampling/pricing on Wednesday (wdate)
% Initiate investment one business day after pricing (wdaten)
% Delta hedge with daily update
% Liuren Wu, liurenwu@gmail.com
%%%%%%%%%%%%
format compact; clearvars;
load('./data/Sampledailydata.mat');
load('./data/Momentsw.mat');

%Expand d1 and maturity index to match week data dimension
d1m=NaN(Tw,nm,nk,nc); matm=d1m;
for m=1:nm
    d1m(:,m,:,1)=repmat(d1(:)',Tw,1);d1m(:,m,:,2)=repmat(d1(:)',Tw,1);
end
for k=1:nk
    matm(:,:,k,1)=repmat(mat(:)',Tw,1);matm(:,:,k,2)=repmat(mat(:)',Tw,1);
end

F0=100; %normalize to 100 forward

NK=NaN(Tw,nm,nk,nc); %normalized strike of initiated contracts
LFK=NK; %lnF/K

PLD=NaN(Tw,nm,nk,nc); %daily delta hedged PL
CPL=PLD; %Call PL
for t= 1:Tw-1
    %Strike, value, delta of contracts at initiation (t)
    for c=1:nc
        for k=1:nk
            xo=Iwn(t,:,k,c);tvo=mat.*xo(:).^2;
            LFK(t,:,k,c)=d1(k)*sqrt(tvo)-0.5*tvo; %ln F/K
        end
    end
    NKt=exp(-LFK(t,:,:,:))*F0;%normalized strike
    NK(t,:,:,:)=NKt;
    
    VT0=Iwn(t,:,:,:).*sqrt(matm(t,:,:,:));
    d10=d1m(t,:,:,:);
    d20=d10-VT0;
    C0=F0.*normcdf(d10)- NKt.*normcdf(d20);
    D0=normcdf(d10);
    
    %Changes over the next week, t+1
    sdate=wdaten(t);
    tn=t+1;
    edate=wdaten(tn);
    tt=find(ddate>sdate&ddate<=edate);
    ndd=length(tt);
    ddt=ddate(tt);
    dLFt=NaN(1,nm,nk,nc,ndd);
    IVKt=dLFt;
    DNt=dLFt;
    CNt=dLFt;
    FNt=dLFt;
    for td=1:ndd
        tn=tt(td);
        dt=(ddt(td)-sdate)/365;
        matdt=mat-dt; %needs to make sure this is greater than 0, true in weekly return
        for c=1:nc
            LFSF=interp1([0;mat], [0;LFS(tn,:,c)'],matdt,'linear'); 
            %LFS in the future for today's contract, with shrinking maturity
            dLFS=LFSF -LFSwn(t,:,c)';
            dlf=log(spot(tn,c)./Swn(t,c))+ dLFS; %dlnF
            
            %interpolate across maturity, linear on total variance
            ivk1=NaN(nm,nk);lfk1=ivk1; lfkp=ivk1;
            for k=1:nk
                dLFt(1,:,k,c,td)=dlf;
                
                %K relative to new F, log(F(t+n)/K=lnF/K +lnF(t+n)/F(t)
                lfkp(:,k)=LFK(t,:,k,c)'+dlf;
                
                %new vol at d1--> interpolate to shrank maturity,
                %interpolation: linear in total variance
                xn=IV(tn,:,k,c);  tvn=mat.*xn(:).^2;
                iim=isfinite(tvn); ni=sum(iim);
                tv1d=NaN(nm,1); %only generate numbers on available maturity
                tv1d(iim)=interp1([0;mat(iim)], [0;tvn(iim)],matdt(iim),'linear');
                ivk1(:,k)=sqrt(tv1d./matdt);
                lfk1(:,k)=d1(k)*sqrt(tv1d)-0.5*tv1d; %relative strike at d1
                
            end
            %interpolate on d1 at each maturity, cubic in Nd1
            for m=1:nm
                ym=ivk1(m,:)'; iik=isfinite(ym); ni=sum(iik); %number of available strikes
                if ni>=3 %at least 3 strikes
                    %extrapolate on Nd1 is more stable because it is within [0,1]
                    ym=ym(iik);
                    tvm=matdt(m)*mean(ym).^2; %tvm is juts used for transformation as a scalar
                    Nd11=normcdf(((lfk1(m,iik)'+.5*tvm)./sqrt(tvm)));
                    
                    Nd1p=normcdf(((lfkp(m,:)'+.5*tvm)./sqrt(tvm)));
                    iiz=isfinite(Nd1p);
                    IVKt(1,m,iiz,c,td)=interp1(Nd11,ym,Nd1p(iiz),'pchip','extrap');
                end
            end
        end
        
        dLFd=dLFt(1,:,:,:,td);
        IVKd=IVKt(1,:,:,:,td);
        
        matdtt=max(0,matm(t,:,:,:)-dt);
        Fn=F0*exp(dLFd);
        VTn=IVKd.*sqrt(matdtt);
        d1n=(log(Fn./NKt)+0.5*VTn.^2)./VTn;
        d2n=d1n-VTn;
        CN=Fn.*normcdf(d1n)- NKt.*normcdf(d2n);
        DN=normcdf(d1n);
        %ii=matdtt==0; %expired options (does not happen here)
        %CN(ii)=max(Fn(ii)-NKt(ii),0);   DN(ii)=sign(CN(ii));
        DNt(1,:,:,:,td)=DN;
        CNt(1,:,:,:,td)=CN;
        FNt(1,:,:,:,td)=Fn;
    end
    
    %aggregate to weekly
    CNN=CNt(:,:,:,:,end);
    DNN=DNt(:,:,:,:,1:end-1);
    
    %hedging PL, sum daily hedging PL
    PLDH=-D0.*(FNt(1,:,:,:,1)-F0);
    for j=2:ndd
        PLDH=PLDH-DNN(:,:,:,:,j-1).*(FNt(1,:,:,:,j)-FNt(1,:,:,:,j-1));
    end
    CPL(t,:,:,:)=(CNN-C0);
    PLD(t,:,:,:)=(CNN-C0)+PLDH;
    
end


save('./data/PLw.mat','wdate','PLD','-mat');


return

